{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Synchronization\n",
    "\n",
    "## Introducation \n",
    "\n",
    "We have discussed how to transmit digitally over the air, utilizing a digital modulation scheme like QPSK and by applying pulse shaping to limit the signal bandwidth. Channel coding can be used to deal with noisy channels, such as when you have low SNR at the receiver. Filtering out as much as possible before digitally processing the signal always helps. In this chapter we will investigate how synchronization is performed on the receiving end. Synchronization is a set of processing that occurs before demodulation and channel decoding. The overall tx-channel-rx chain is shown below, with the blocks discussed in this chapter highlighted in yellow. (This diagram is not all-encompassing–most systems also include equalization and multiplexing).\n",
    "\n",
    "![Alt text](https://pysdr.org/_images/sync-diagram.svg)\n",
    "\n",
    "## Simulating Wireless Channel\n",
    "\n",
    "Before we learn how to implement time and frequency synchronization, we need to make our simulated signals more realistic. Without adding some random time delay, the act of synchronizing in time is trivial. In fact, you only need to take into account the sample delay of any filters you use. We also want to simulate a frequency offset because, as we will discuss, oscillators are not perfect; there will always be some offset between the transmitter and receiver’s center frequency.\n",
    "\n",
    "The Python code in this chapter will start from the code we wrote during the pulse shaping Python exercise: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import signal\n",
    "import math\n",
    "\n",
    "# this part came from pulse shaping exercise\n",
    "num_symbols = 100\n",
    "sps = 8\n",
    "bits = np.random.randint(0, 2, num_symbols) # Our data to be transmitted, 1's and 0's\n",
    "pulse_train = np.array([])\n",
    "for bit in bits:\n",
    "    pulse = np.zeros(sps)\n",
    "    pulse[0] = bit*2-1 # set the first value to either a 1 or -1\n",
    "    pulse_train = np.concatenate((pulse_train, pulse)) # add the 8 samples to the signal\n",
    "\n",
    "# Create our raised-cosine filter\n",
    "num_taps = 101\n",
    "beta = 0.35\n",
    "Ts = sps # Assume sample rate is 1 Hz, so sample period is 1, so *symbol* period is 8\n",
    "t = np.arange(-51, 52) # remember it's not inclusive of final number\n",
    "h = np.sinc(t/Ts) * np.cos(np.pi*beta*t/Ts) / (1 - (2*beta*t/Ts)**2)\n",
    "\n",
    "# Filter our signal, in order to apply the pulse shaping\n",
    "samples = np.convolve(pulse_train, h)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Adding a Delay\n",
    "We can easily simulate a delay by shifting samples, but it only simulates a delay that is an integer multiple of our sample period. In the real world the delay will be some fraction of a sample period. We can simulate the delay of a fraction of a sample by making a “fractional delay” filter, which passes all frequencies but delays the samples by some amount that isn’t limited to the sample interval. You can think of it as an all-pass filter that applies the same phase shift to all frequencies. (Recall that a time delay and phase shift are equivalent.) The Python code to create this filter is shown below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create and apply a fractional delay filter\n",
    "delay = 0.4 # Fractional delay in samples\n",
    "N = 21 # Number of taps\n",
    "n = np.arange(-N//2, N//2)\n",
    "plt.plot(np.sinc(n))\n",
    "h = np.sinc(n - delay)\n",
    "plt.plot(h)\n",
    "plt.show()\n",
    "h *= np.hamming(N)\n",
    "plt.figure()\n",
    "plt.plot(h)\n",
    "h /= np.sum(h) # Normalize to avoid chaning amplitude\n",
    "plt.figure()\n",
    "plt.plot(h)\n",
    "samples = np.convolve(samples, h)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, we are calculating the filter taps using a sinc() function. A sinc in the time domain is a rectangle in the frequency domain, and our rectangle for this filter spans the entire frequency range of our signal. This filter does not reshape the signal, it just delays it in time. In our example we are delaying by 0.4 of a sample. Keep in mind that applying any filter, delays a signal by half of the filter taps minus one, due to the act of convolving the signal through the filter."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Adding a Frequency Offset\n",
    "\n",
    "To make our simulated signal more realistic, we will apply a frequency offset. Let’s say that our sample rate in this simulation is 1 MHz (it doesn’t actually matter what it is, but you’ll see why it makes it easier to choose a number). If we want to simulate a frequency offset of 13 kHz (some arbitrary number), we can do it via the following code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# apply a freq offset\n",
    "fs = 1e6 # assume our sample rate is 1 MHz\n",
    "fo = 13000 # simulate freq offset\n",
    "Ts = 1/fs # calc sample period\n",
    "t = np.arange(0, Ts*len(samples), Ts) # create time vector\n",
    "samples = samples * np.exp(1j*2*np.pi*fo*t) # perform freq shift"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Time Synchronization"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following Python code implements the Mueller and Muller clock recovery technique."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu = 0\n",
    "out = np.zeros(len(samples) + 10, dtype=complex)\n",
    "# stores values, each iteration we need the \n",
    "# previous 2 values plus current value\n",
    "out_rail = np.zeros(len(samples) + 10, dtype=complex)\n",
    "i_in = 0 # input samples index\n",
    "i_out = 2 # output index (let first two outputs be 0)\n",
    "while i_out < len(samples) and i_in+16 < len(samples):\n",
    "    out[i_out] = samples[i_in + int(mu)] # grab what we think is the \"best\" sample\n",
    "    out_rail[i_out] = int(np.real(out[i_out]) > 0) + 1j*int(np.imag(out[i_out]) > 0)\n",
    "    x = (out_rail[i_out] - out_rail[i_out-2]) * np.conj(out[i_out-1])\n",
    "    y = (out[i_out] - out[i_out-2]) * np.conj(out_rail[i_out-1])\n",
    "    mm_val = np.real(y - x)\n",
    "    mu += sps + 0.3*mm_val\n",
    "    i_in += int(np.floor(mu)) # round down to nearest int since we are using it as an index\n",
    "    mu = mu - np.floor(mu) # remove the integer part of mu\n",
    "    i_out += 1 # increment output index\n",
    "out = out[2:i_out] # remove the first two, and anything after i_out (that was never filled out)\n",
    "plt.plot(out)\n",
    "plt.show()\n",
    "samples = out # only include this line if you want to connect this code snippet with the Costas Loop later on\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.4"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
